%
% Tip of the Iceberg : Generate simulated dataset with measurement error, and then estimate model
%
% moxnes@gmail.com. August 2014.

%
% (1) Given a matrix of random draws, solve model, and construct simulated dataset
%

% Parameters
clear;

workdir = '/Dropbox/RnD_KK/tradecost_replication';
% addpath1 = strcat(workdir,'/matlab_model2/t_tau2');
addpath2 = strcat(workdir,'/ReStat_code/func');

% path(addpath1,path);
path(addpath2,path);

params.K = 4;  % Countries
K = params.K;
params.gdp = ones(params.K,1);

% Standard deviation of demand shock
d_stdev = .92;  

% Number of bootstraps
bootn = 150;
sigma = [3 4]';
params.gamma = 4;                 % Simonovska Waugh

params.phi = 0;
params.P = 2;                      % Number of products
params.draws = 5000;               % No. of random draws         

% Country size
params.L = params.gdp;    
params.mu = 1;                    % Consumption share for monopolistic good

% Random productivity draws. Productivity x is Pareto, F(x) = 1-x^(-gamma)
% Bounded from [1,+inf>, costs are <0,1]

rand('seed',1);
r = rand(params.draws,1);
params.x = (1-r).^(-1/params.gamma);

% Fixed costs product #1 (high sigma)
cxforeign = 0.94;  
cxhome = 0; 
Cxall(:,:,1) = ones(params.K)*cxforeign-eye(params.K)*cxforeign+eye(params.K)*cxhome;   % (KxK)

% Fixed costs product #1 (low sigma)
cxforeign = 0.35;  % 
Cxall(:,:,2) = ones(params.K)*cxforeign-eye(params.K)*cxforeign+eye(params.K)*cxhome;   % (KxK)

kappa = .12; 
params.potential_entrants = params.L*kappa;   
% set this so that home cutoff=1

% Product specific trade cost
t_p = [.8 1]';

% Destination specific trade cost
t_n = [0 .1 .2 .3];
 
disp('lnt_n'); disp(log(t_n(2:end)));
disp('lnt_p'); disp(log(t_p'));

% --------------------------------------
% Solve the model X times, with different random draws each time
% --------------------------------------
for vv = 1:bootn
disp('Bootstrap #'); disp(vv);
% Lognormal demand shocks, per market
%rand('seed',2);
params.d = randn(params.draws,K);

% ------------------------------------------------------
for pp=1:length(t_p)  % Loop over products
    
% Additive barriers
tmp = repmat(t_n,K,1)*t_p(pp);
params.t = tmp - diag(diag(tmp));

% Multiplicative barriers matrix, cols=to, rows=from
const = 1.1;
params.Tau = const*ones(K) - (const-1)*diag(diag(ones(K)));

params.sigma = sigma(pp);
s = params.sigma;
markup = s/(s-1);
params.delta1 = (s/params.mu)^(1/(s-1))*(s/(s-1));

% Fixed cost for this product
params.Cx = Cxall(:,:,pp);

% ---------------------------------
% Check parameters before we move on
% ---------------------------------

%if (params.gamma<params.sigma-1) error('Error: sigma too high'); end;
if (~all(all(params.Tau.^(params.sigma-1).*params.Cx >= repmat(diag(params.Cx),1,params.K))))
    error('Error: Too low export costs');
end;

% What is the price index if x = mean(x)?
% Use the P0 as starting value.
p = markup.*(params.Tau./mean(params.x)+params.t);

for n=1:params.K  
  P0(n,1) = sum(params.potential_entrants.*p(:,n).^(1-s)).^(1/(1-s));
end

%P0 = P0*2;      % Start with a higher P0
if (any(isnan(imo_solve(P0,params)))) error('Cannot solve price index. Problem ill conditioned.'); end;

% -----------------------------
% Solve equilibrium
% -----------------------------

func = @(p) imo_solve(p,params);

% % Solve f(x)=0
options = optimset('Display','none','MaxIter',1000,'MaxFunEvals',20000);
tic; [P1,fval,exitflag,output] = fsolve(func,P0,options); toc;
if exitflag<=0 error('fsolve does not converge'); end;

% Find x_x and pi simulaneously
pi = solve_pi_xx(P1,params);
if (pi==0) error('Profits negative, set to zero'); end;

Y = params.L*(1+pi);                     % y is total income
x_x = cutoff_imo_x(P1,Y,params);
          
    i=1;  % Choose country 1 as source
    for n=1:K
      export = params.x>x_x(i,n);
      firms(n,pp) = sum(export);
      if firms<200 warning('Less than 200 firms exporting'); end;      
      pricetmp = markup.*(params.Tau(i,n)./params.x+params.t(i,n));
      pr(:,n,pp) = export.*pricetmp;
      prtilde(:,n,pp) = pr(:,n,pp)/params.Tau(i,n) - params.t(i,n)/params.Tau(i,n); 
      dshock = exp(d_stdev*params.d(:,n));
      qty(:,n,pp) = export.*params.mu*Y(n).*pricetmp.^((s-1)*params.phi-s).*dshock.*P1(n).^(s-1);
      sales(:,n,pp) = qty(:,n,pp).*prtilde(:,n,pp);      
    end
    
    % CALIBRATE Share exporting firms (for home country=1)
    homef = params.x>x_x(1,1);
    exportf = any(repmat(params.x,1,K-1)>repmat(x_x(1,2:end),params.draws,1),2); 
    ex_share = sum(exportf)/sum(homef);
    ex_share_all(vv,pp) = ex_share;
    %disp('Share exporting firms'); disp(ex_share);    
end
    
%
% Create simulated dataset. Export firm-level quantity and price export data. Source=1
%
    
% file1: firmID productID qty1 qty2.. qtyK (sorted by productID)
% file2: firmID productID fobprice1 fobprice2 .. fobpriceK
% file3: productID number_of_firms


firmarr = (1:params.draws)';
prodarr1 = ones(params.draws,1)*1;
prodarr2 = ones(params.draws,1)*2;                         
       
% ----------------------------------
% Output from model
% ----------------------------------

% Trade costs relative to fob price
for n=1:K
  tmp = prtilde(:,n,1);
  prtilde_median(n,1) = median(tmp(tmp>0));
  tmp = prtilde(:,n,2);
  prtilde_median(n,2) = median(tmp(tmp>0));
end

tmat = repmat(t_n',1,params.P).*repmat(t_p',params.K,1);
taumat = repmat(params.Tau(1,:)',1,params.P);

trueTC = (tmat./taumat)./prtilde_median;
disp('MODEL Trade costs t/tau relative to median export price');
disp(trueTC);

% --------------------------------------------------
% Add measurement error to prices and quantities
% --------------------------------------------------

stderr = std(prtilde(:))*0.1;


% Standard deviations of measurement error
rand('seed',10); err1 = randn(params.draws,params.K,2)*.06;

stderr
std(exp(err1(:))) % Should be similar to stderr

% prtilde = prtilde.*exp(err1);

prtildeEX = [firmarr prodarr1 prtilde(:,2:end,1);
             firmarr prodarr2 prtilde(:,2:end,2)];
         
qtyEX = [firmarr prodarr1 qty(:,2:end,1); 
         firmarr prodarr2 qty(:,2:end,2)];

countEX = [1 params.draws; 
           2 params.draws];  


% --------------------------------------------------- 
% (2) Given simulated data, estimate model 
% --------------------------------------------------- 

% Load data
M = qtyEX;
M2 = prtildeEX;
data.numfirms = countEX;

thresh = '40';

data.P = size(data.numfirms,1);
data.K = size(M,2)-2;           % # countries
thresh2 = str2num(thresh);
start = 1;
for p=1:data.P
  stop = start+data.numfirms(p,2)-1;
  qty = M(start:stop,3:end);    % I_k x K Sales (quantity)
  price = M2(start:stop,3:end);
  entrants  = qty>0;            % (I_k x P) x K Export participation
  
  for n=1:data.K
    tmp = qty(:,n);
    tmp_p = price(:,n);
    if length(tmp(entrants(:,n)))>=thresh2           % Check if we want to use this product-dest pair
      data.usedata(n,p)=1;
      data.sales_actual{n,p} = tmp(entrants(:,n));
      data.lnX{n,p} = log(data.sales_actual{n,p});
      data.price{n,p} = tmp_p(entrants(:,n));     
      data.price_mean(n,p) = mean(data.price{n,p});
      data.price_median(n,p) = median(data.price{n,p});
      data.entrants(n,p) = length(data.sales_actual{n,p});      
    else data.usedata(n,p)=0; end;
  end  
  start = stop+1;
end

% Number of destinations per product
data.K2 = sum(data.usedata)'; 
data.avgentrants = mean(data.entrants(logical(data.usedata)));
disp('Average # entrants, across destinations'); disp(data.avgentrants); 

clear M N;
res = []; lb = []; ub = [];

res = [];
options = optimset('Display','none','TolFun',1e-8,'TolX',1e-8);

    for p=1:data.P
          prm.p = p;
          disp('Product #'); disp(p);
          disp('Number of destinations'); disp(data.K2(p));
          f = @(b) imo_2011_NLSFE(b,data,prm,0);
          
          tmp = data.price_median(:,p);
          tmp = tmp(tmp>0);     % Median prices from destinations with positive sales
          nn = length(tmp);
          b0 = [log(tmp'*.1) ones(1,nn)*5 -2]; 
          lb = [log(tmp'*.01) ones(1,nn)*-10 -10];          
          ub = [log(tmp'*1) ones(1,nn)*30 10]; 
          % t max 1x price
          
          [b,fval,exitflag,output,lambda,grad,hessian]  = fmincon(f,b0,[],[],[],[],lb,ub,[],options);
          if exitflag<0 error('No estimate'); end
          lnt = b(1:data.K2(p));
          a = b(data.K2(p)+1:2*data.K2(p));
          phisig1 = b(end);

          % Find the correct indices for t
          i = 1; t_arr = []; a_arr = [];         
          for n=1:data.K
            if (data.usedata(n,p)==1) 
                lnt_arr(n) = lnt(i); 
                a_arr(n) = a(i);
                i = i+1; 
            else
                lnt_arr(n) = NaN;
                a_arr(n) = NaN;

            end
          end
          res.lntmat(:,p) = lnt_arr';
          res.amat(:,p) = a_arr';
          res.phisig1(p) = phisig1;
    end
        
    % Collect results in a matrix
    phisig1_pmat(vv,:) = res.phisig1;
    lntmatall(:,:,vv) = res.lntmat;
    tc_mat(:,:,vv) = exp(res.lntmat)./data.price_median;    
end


disp('--------- DONE ------------');

%%% ........................%%%
%%% TABLE 4                 %%%
%%% ........................%%%

disp('TRUE: TC'); disp(trueTC); 
disp('TRUE: sigma'); disp(sigma'); 


tc_mat2 = exp(mean(lntmatall,3))./data.price_median;
disp('ESTIMATES: MEAN TC ALT'); 
disp(tc_mat2);

disp('ESTIMATES: Standard deviation TC');
disp(std(tc_mat,0,3));

disp('phisig1_p mean'); disp(mean(phisig1_pmat));
disp('phisig1_p stdev'); disp(std(phisig1_pmat));

% Entrants
disp('Average # entrants');
kk = firms(2:end,:);
disp(mean(kk(:)));

      
      
%%% ........................%%%
%%% TABLE 5                 %%%
%%% ........................%%%

disp('ESTIMATES: Mean TC'); 
disp(mean(tc_mat,3));

disp('ESTIMATES: Standard deviation TC');
disp(std(tc_mat,0,3));

disp('phisig1_p mean'); disp(mean(phisig1_pmat));
disp('phisig1_p stdev'); disp(std(phisig1_pmat));

disp('Average # entrants');
kk = firms(2:end,:);
disp(mean(kk(:)));
      







% CALIBRATE sddev ln eps relative to mean sales 
for kk=1:2
    for nn=1:3 
        dat = data.sales_actual{nn,kk};
        epsRel(nn,kk) = d_stdev/mean(log(data.sales_actual{nn,kk}.*data.price{nn,kk}));
        stdevSales(nn,kk) = std(log(data.sales_actual{nn,kk}));
        stdevSales2(nn,kk) = std(log(data.sales_actual{nn,kk}.*data.price{nn,kk}));
        lnpratio(nn,kk) = log(prctile(dat,90)/prctile(dat,10));
    end
end
disp('stddev log qty'); disp(stdevSales);
disp('log p90/10 ratio'); disp(lnpratio);
disp('median'); disp(median(lnpratio(:)));

% CALIBRATE mean share exporting firms
disp('Mean share exporting firms');
disp(mean(ex_share_all));
disp(mean(ex_share_all(:)));

% Entrants
disp('Average # entrants');
kk = firms(2:end,:);
disp(mean(kk(:)));
      